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We report numerical calculations of the concentration of interstitials in hard-sphere crystals. We 
find that, in a three-dimensional fee hard-sphere crystal at the melting point, the concentration of 
interstitials is 2.7(4) • 1CP 8 . This is some three orders of magnitude lower than the concentration of 
vacancies. A simple, analytical estimate yields a value that is in fair agreement with the numerical 
results. 



I. INTRODUCTION 

Any crystal in equilibrium will contain defects, such as vacancies, interstitials and dislocations. Of these, the point 
defects are the most common. Some thirty years ago, Bennett and Alder jlj estimated the equilibrium concentration 
of vacancies in a hard-sphere crystal and found that, close to melting, this concentration could be quite high (typically, 
one vacancy per 4000 lattice sites). At present, the question of the concentration (and transport) of point defects 
in (colloidal) hard-sphere crystals takes on a renewed — and now quite practical — significance. Apart from the 
theoretical interest in hard sphere crystals as a model system, crystals from colloidal particles, having lattice sizes 
comparable to the wavelength of light, are being prepared and studied because of their potentially interesting photonic 
properties. Clearly the presence of even a small number of defects can have a pronounced effect on the nature of 
photonic states in such materials. Moreover, as the accuracy of free energy calculations increases, it is no longer 
permissible to ignore the contribution of vacancies to the total free energy. The aim of the present paper is to review 
briefly the statistical mechanical description of a crystal with point defects. This problem is not completely trivial, 
as the concept of a vacancy or interstitial is inextricably linked to that of lattice sites. And lattice sites loose their 
meaning in a disordered state. So, we should first address the question: when is it permissible to count states with a 
different number of lattice sites as distinct? The answer is, of course, that this is only true if these different states can 
be assigned to distinct volumes in phase space. This is possible if we impose that every particle in a crystal is confined 
to its Wigner-Seitz cell. In three dimensional crystals, this constraint on the positions of all particles has little effect 
on the free energy (in contrast, in a liquid it is not at all permissible). In a two-dimensional crystal, the constraint 
is more problematic, at least in the thermodynamic limit. However, for large but finite two-dimensional systems, the 
single-occupancy cell constraint is also quite reasonable. Below, we describe two alternative (but equivalent) routes to 
arrive at the free energy of a crystal with vacancies. In one case, we use the Grand-Canonical ensemble. This would 
seems to be the most obvious ensemble to use when describing a system with a fluctuating number of particles. Yet, 
the analysis is complicated by the fact that not only the number of particles, but also the number of lattice sites, 
may fluctuate. In the second analysis, we consider an isothcrmal-isobaric system. The latter approach is simpler and 
is, apart from a minor correction, equivalent to the one followed by Bennett and Alder jjj. We then describe our 
numerical approach to compute the concentration of interstitials in a hard-sphere crystal. We compare our numerical 
results with a simple theoretical estimate. 

II. FREE ENERGY OF VACANCIES 

A. The grand-canonical route 

When considering the statistical mechanics of a crystal with vacancies, it is convenient to consider first a system 
with a fixed number of lattice sites, M , contained in a volume V . If this crystal is in contact with a particle-reservoir 
at chemical potential /i, then the number of vacancies in the crystal may fluctuate. In principle, the crystal could 
also contain interstitials but, for the time being, we shall ignore this possibility. It is then easy to write down the 
expression for the grand potential of the crystal E': 
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where /3 = l/k B T. Note that this is not the true grand potential, because we should also allow for fluctuations in the 
number of lattice sites. We denote the free energy of a crystal with no vacancies by F<°) = -kT\nQ M - In practice, 
the equilibrium concentration of vacancies in a crystal is very low. We shall therefore make the approximation that 
vacancies do not interact. This assumption is not as reasonable as it seems, as the interaction of vacancies through 
the stress field is quite long-ranged. The assumption that vacancies are ideal makes it easier to compute the canonical 
partition function of a crystal with n vacancies: 
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where we have used the notation f( n ) to denote the free energy of a crystal with n vacancies at given positions. As 
the vacancies are assumed to be non-interacting, it is clear that we can write 



F (n) = p(0) _ nfi = Mfo _ nfi 



(3) 



where fo is the free energy per particle in the defect-free crystal and -/i is the change in free energy of a crystal due 
to the creation of a single vacancy at a specific lattice point[]. Combining Eqns. ([j]), (||) and (H), we obtain 
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Usually, exp[— /3(/i — /i)] is much less than unity. This allows us to write: 
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the average number of vacancies follows as 
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Now, we should take into account the fact that, actually, the number of lattice sites itself is not fixed but will adjust 
to the number of vacancies. The total grand partition function is therefore a sum over all states with different number 
of lattice sites M' = M + AM. In practice, AM <§; M. We can then write 

oo 

AA/=-oo 

Note that 5' depends on both M + AM and \i. We now choose the reference number of lattice sites such that in that 
particular case y, is equal to the chemical potential of the perfect lattice. That is: 

F (o) +P (o) v = M/i 

We also introduce a rather strange quantity, namely the "grand canonical" partition function of a perfect lattice with 
a fixed number of particles (M) : = Qm exp(M/3/i). Of course, 3^ is not a true grand-canonical partition 
function, as the number of particles in this system is fixed. The Grand Canonical partition function then becomes: 



1 The choice of the minus sign in the definition will later turn out to be convenient. 
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We assume (as usual) that, in the thermodynamic limit, 5 is dominated by the largest term in the sum. Hence, we 
have to determine the point where the derivative of 5 with respect to M vanishes. To this end, we perform a Taylor 
expansion of the exponent in Eqn. (j7j) in powers of AM. Note that for a perfect — defect free — lattice 
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By combining Eqn. (§)-© we obtain 
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Note that /i = F<® - P (1) . Hence, 
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where AP 0,1 is the difference in the pressure of two crystals with M lattice sites, one with zero vacancies and the 
other with one vacancy (both in a fixed volume V and a temperature T). As a consequence, 



(M + AM)exp - P([i - ft [M + AM] ) 

« Mexp { - (3[p ~ /i(M)]}|l + + (/3AP -Vp)Am| 
Inserting the Taylor expansion in the expression for S, we obtain 
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It may seem inconsistent that we expand to second order in AM in the last term of Eqn. ( |l4|) but only to first order 
in the preceding terms. However, as we show below, we actually expand — consistently — to second order in the 
vacancy concentration. 

We define the fractional change in the number of lattice sites, y, as 
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Similarly, we define the vacancy concentration x as 
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Finding the maximum in the exponent is then equivalent to maximizing 
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where x(0) is the value of x for AM — and p s is the density of the ideal reference lattice. The maximum value is at 
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and the maximum of the grand canonical potential is 
S ps exp(-/3F (0) (M) + Mpjj) exp 
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Hence, the presence of vacancies increases the grand potential (as it should) and it changes (increases) the number of 
lattice sites. However, ignoring terms of order 0(x 2 (O)), the vacancy concentration is still given by Eqn. (||). 

The next question is: how do vacancies affect the melting curve. Now our definition of the reference system (i.e. 
the perfect lattice with the same fi) turns out to be convenient. Note that the Grand Canonical partition function is 
related to the pressure by 

I3PV = lnS = (3P (0) V + Mx(0) +O{x(0) 2 ) (20) 

Ignoring terms quadratic in the vacancy concentration, we find that the effect of allowing for vacancies is to increase 
the pressure of the solid by an amount 

AP w x(0)p s kT (21) 

Let us assume that the liquid was in equilibrium with the perfect crystal at pressure P and chemical potential /i. 
Then it is easy to verify that the shift in the coexistence pressure due to the presence of vacancies is 

_ -x(0)kT 

O^coex _i _i V / 

Pi - Ps 

and the corresponding shift in the chemical potential at coexistence is 

OfJ-coex — \^"J 
Pi 

Direct calculations of the vacancy concentration in a hard-sphere crystal at melting [Q indicate that a;(0) pa 2.6 • 10~ 4 . 
Hence, the increase in the coexistence pressure due to vacancies is SP coex pa — 2.57- 10~ 3 kT/a 3 (where a is the particle 
diameter). The corresponding shift in the chemical potential at coexistence is Sp, coex = —2.74- 10~ 3 A:T Note that this 
shift is very significant when compared to the accuracy of absolute free-energy calculations of the crystalline solid [|| . 



B. The NPT route 

Bennett and Alder |l| work with the Gibbs free energy rather than the Helmholtz free energy. Their expression for 
the vacancy concentration is based on the analysis of the effect of vacancies on the Gibbs free energy of a system of 
N particles at constant pressure and temperature. First, we define g vac , the variation in the Gibbs free energy of a 
crystal of M particles due to the introduction of a single vacancy at a given lattice position 

g vac = G M+ i,i(N,P,T) - G M>0 (N,P,T) 

— Fm+i,i(Vm+i,i) — -Fm,o(Vm,o) + P{Vm+i,i — Vm,o) (24) 

where the first subscript refers to the number of lattice sites in the system, and the second subscript to the number 
of vacancies. In this equation we distinguish M, the original number of lattice sites, and N , the number of particles, 
even though in the present case N = M. Let us write j\ (Eqn. (|37|)) as 

—fx = Fm+ia(Vm+i,o) — Fm+i,o(Vm+i,o) 
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Hence, 

gvac _ Fm+1,i(Vm+1,i) — ^M+1,i(Vm+1,o) + 

Fm+i,i{Vm+i,o) — Fm+i,o(Vm+i,o) + 
^M+i,o(Vjf+i,o) — Fm,o(Vm,o) + 

P(V m +ia - V Mt0 ) (25) 

The next step is to introduce a hypothetical defect free crystal with M lattice sites, at the same pressure as the system 
with M + 1 lattice sites. The volume of this system is Vm,o — {M/(M + l)}Vjf + i ) o. Similarly, the free energy of this 
hypothetical system is Fm,o = {M/(M + 1)}Fm+i,o- Note also that 

F m +i.i{Vm+i,i) - F M ,o{V M fi) = -PAv - h + fo (26) 

where Av = v vac — v part is the difference in volume of a vacancy and a particle, at constant pressure and number of 
lattice sites. Moreover, 

P(V M+ i.i - V M ,o) = P(Av + V/N) (27) 

Hence, the Gibbs free energy difference associated with the formation of a vacancy at a specific lattice site, Gm,i — 
Gm-i,o = 9 vac i is then 

9 vac = P(V M+ i.i - Vm.o) - h + (Av + V/N)P + f 
= P{Vm,i — Vm,o + Vm,o — Vm-i,o) — fi + fo 
= P(V/N) - h + fo 
= {P/p + fo)-h 

= Mo - h (28) 

where we have defined no = (P/p + fo)- Now we have to include the entropic contribution due to the distribution of 
n vacancies over M lattice sites. This total Gibbs free energy then becomes 

G = G a (N) + ng ^ + MkT In ± + [l - £] In [l - ^] ) (29) 

n 

w Go (N) + ng vac + nkT In — - nkT (30) 
If we minimize the Gibbs free energy with respect to n, we find 

(n) McM~P9 vac ) 

where we have ignored a small correction due to the variation of InM with n. If we insert this value in the expression 
for the total Gibbs free energy, we find: 

G = G (AO + (n) g vac - (n) g vac - (n) kT = G - (n) kT 

The total number of particles is M— (n). Hence the Gibbs free energy per particle is 

Go - (n) kT (n) kT 

/i= N =^°-^ 

w mo - x v kT (31) 

Thus the change in chemical potential of the solid is 

Ap = -x v kT (32) 
from which it follows that the change in pressure of the solid at fixed chemical potential is equal to 

AP = x vPs kT (33) 
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This is equivalent to Eqn. (Ell) above. Hence the Bennett- Alder scheme is equivalent to the "grand-canonical" scheme 

It should be pointed out that the variation in volume due to the replacement of a particle by a vacancy can be 
computed either directly, in a constant-pressure simulation, or indirectly by measuring the change in pressure in a 
constant volume simulation. The two methods are related through the thermodynamic relation 



dV\ { dP\ {dN 



dP J NT \ dN J VT \ dV / PT 



= -1 (34) 



Noting that the number of vacancies is n = M — N, we can see that the change in pressure with the number of 
vacancies for a fixed number of lattice sites, is 

dP\ {dP\ fdV 



dN) VT \dVj NT \dN /pr 



(35) 

In particular, the pressure change due to one vacancy (i.e. AP ' 1 , defined in Eqn. (Oh) is 



III. COMPUTATIONAL SCHEME 
A. Vacancies 

Numerically, it is straightforward to compute the equilibrium vacancy concentration. As before, the central quantity 
that needs to be computed is — fx, the change in free energy of a crystal due to the creation of a single vacancy at 
a specific lattice point. In fact, it is more convenient to consider +/i, the change in free energy due to the removal 
of a vacancy at a specific lattice point. This quantity can be computed in several ways. For instance, we could use a 
particle-insertion method. We start with a crystal containing one single vacancy and attempt a trial insertion in the 
Wigner-Seitz cell surrounding that vacancy. Then f\ is given by 

/l = - ff h(^l^) (37) 

where Vws 1S the volume of the Wigner-Seitz cell, and A 17 is the change in potential energy associated with the 
insertion of a trial particle. For hard particles 

A = -kTln( ^ WsPa ™} Vws) } (38) 

where P aC c(Vws) is the probability that the trial insertion in the Wigner-Seitz cell will be accepted. As most of the 
Wigner-Seitz cell is not accessible, it is more efficient to attempt insertion in a sub-volume (typically of the order of 
the cell- volume in a lattice-gas model of the solid). However, then wc also should consider the reverse move — the 
removal of a particle from a sub- volume v of the Wigner-Seitz cell, in a crystal without vacancies. The only thing 
we need to compute in this case is P rem (v), the probability that a particle happens to be inside this volume. The 
expression for /i is then 

Of course, in the final expression for the vacancy concentration, the factor A d drops out (as it should), because it is 
cancelled by the same term in the ideal part of the chemical potential. 



In Ref. M, a slightly different expression is found, but this is due to a small error in the derivation in that paper. 



G 



B. Interstitials 



As in the case of vacancies, the calculation of interstitials centers around the calculation of <?/, the free energy 
associated with introducing an interstitial into the system (in the NPT ensemble), gi can be expressed as a sum of 
two parts: the free energy of introducing a point particle into the system (g ms — — fcTln < 1 — rj >, where rj is the 
packing fraction), and the free energy of growing that particle to the same diameter (erg = 1) as the other particles 
(gB r0W \ 

For the calculation of g9 row , we simulate an extended system consisting of a lattice of N particles with diameter 
(To = 1 on (or near) lattice sites and one extra (interstitial) particle that has a diameter 07 that can vary freely. We 
interpret 07 as an additional coordinate and, in this sense, the system that we are considering is an extended system. 
The partition sum for the extended system is 



Q(N+ l,P,T) = [ da' Q(N + l,P, T, a') 
Jo 



(40) 



where Q(N + 1, P, T, a') is the partition function for the isothermal-isobaric system with one interstitial particle with 
radius a'. The probability of finding the interstitial particle with a specific radius 0' = 07 is 

PioAN + 1 P T)= S 1 ^Q{N + l,P,Ty)8(a I -o') Q(N + 1, P, T, *j) 
( 11 1 ' ' Q(N,P,T) Q(N + 1,P,T) 1 j 

and the Gibbs free energy G(N + 1, P, T, 07) of a system with an interstitial with diameter 07 is equal to 

G(N + 1, P, T, aj) = -k B T\nQ(N + 1, P, T, aj) (42) 

Thus, the free energy difference between a system with a pointlike interstitial (07 = 0) and a full-grown interstitial 
(a-/ = 1) is 

g grow = Q( N + ^ p T _ ^ _ Q( N + Xj p p^ 

Q(N + 1,P,T,1) 



-fcTln 



Q(N + l,P,T,0) 



P(0\N+1,P,T) 

= kTln P (l\N+l,P,T) (43) 

It is obvious that P(a\N + 1,P, T) is can be very small for large values of 07. In order to get an accurate 
histogram for P(<tj\N + 1,P, T), we have to use a biased sampling scheme. We employ the method of umbrella- 
sampling/multicanonical sampling where we associate a weight £(ct) with a, which we use while sampling over 

a: 

P{a\N + 1, P, T, {£}) cx P(a\N + 1, P, T)e^ a) 

If we sample over this distribution, we get a histogram P({a}\N + 1,P, T, {£}), for which we can get the desired 
histogram P({<j}\N + 1, P, T) by refolding the bias: 



P(a\N + 1, P, T) cx P(a\N + 1, P, T, {£}) 



The weights £(cr) are obtained by iteratively running the system and calculating (for the run from the results 

of run i) 

6+i (a) = - In P(<T\N,P,T,{ti(<T)}) +C 

where C is an arbitrary constant |Q. This will make the histogram P({ao}\N + 1,P, T, {£}) converge to a fiat 
distribution over the accessible range. 
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C. Interstitial type discrimination 



The fee crystal has two types of possible places, or 'holes' in which interstitials can reside: one of octahedral shape 
and one of tetrahedral shape. There are four octahedral holes and eight tetrahedral holes in one fee unit cell. To 
measure the relative concentrations of interstitials in these two types of holes, it turns out it is not possible to try to 
prepare the system in one hole and calculate gS row ; because during the course of a simulation the interstitial makes 
many hops. These hops are caused by the interstitial taking the place of a lattice particle which then becomes the 
interstitial. 

In order to measure the relative occupation probability of the different holes, the interstitial has to be traced. This 
is done by the following scheme: at the start of the simulation, every particle i, except the original interstitial, is 
assigned to a lattice position R,. At fixed sampling intervals, the squared distance between the original interstitial 
and the nearest lattice sites (§f nti — {Tint — R-i) 2 ) is compared with <5f = (r* — R^) 2 . If S^ nt i < <5 2 , the interstitial and 
particle i exchange identity (i.e. the interstitial acquires a lattice position Rj and particle i becomes the interstitial). 
Once we have identified the interstitial, it is straightforward to assign it to a tetrahedral or octahedral hole. 

IV. SIMULATION RESULTS 

The free energy calculations were performed with 256 + 1 particle systems (4x4x4 cubic fee unit cells) at four 
different pressures. Different parts of the histogram P({a}\N, P, T, {77}) were calculated in parallel, and subsequently 
combined. The calculation of the weights took about 20 iterations of 2 • 10 4 MC sweeps per CPU on 5 CPU's. Once 
the weights were known for one pressure, they could be used as starting points for the other pressures, accelerating 
the weight calculation considerably. Final calculations were done with approximately 20 iterations of 4 • 10 5 sweeps 
each (again on 5 CPU's). The final P(a\N, P, T) histograms for all four pressures are plotted in Fig. |l|. 

For the calculation of /j at the different pressures, the results for the free energy of the perfect crystal j|] were 
used together with Hall's J|] equation of state. The results are summarized in Fig. ^| and table H. For one pressure 
(P = 11.7, the coexistence pressure), we calculated gi for a larger system (N = 8x8x8 + 1 = 2048 + 1) to check 
for finite-size effects; as can be seen from the results, these are negligible. Using the interstitial type discrimination 
algorithm described above, it was found that the (bigger) octahedral holes are far more likely to contain the interstitial 
than the tetrahedral holes (see table .|). 

V. ANALYTICAL ESTIMATE OF THE FREE ENERGY OF INTERSTITIALS 

As octahedral holes are the largest cavities in a fee crystal, we limit our analysis to these. The number of octahedral 
holes in a fee crystal is equal to the number of lattice sites, the derivation of the expression for the concentration 
of interstitials is almost identical to the one for the vacancy concentration. Let us denote the change in free energy 
associated with the introduction of an interstitial at a specific octahedral site by /j. The concentration of interstitials 
is then 

ar/=earp -/?(//- fi) (44) 

In a static lattice, ro, the radius of such octahedral holes equals (^2/2 — 0.5)0, where a is the nearest-neighbor distance. 
For a hard-sphere crystal at melting, this radius equals ro = 0.229 a. Clearly, in order to fit in an interstitial, the 
cavity has to be expanded. If we assume that the crystal is elastically isotropic (a fair approximation for a cubic 
crystal) then the work need to create a cavity of radius r equals || 

W = 87r/jr (r-r ) 2 (45) 

where fi is the shear Lame coefficient. How large should r be? Clearly, it should be at least 0.5<r, otherwise the 
interstitial would not fit into the lattice. But, in fact, it should be larger, because the interstitial particle itself 
requires some free volume Up- We should therefore minimize the sum of the free energy of a particle in a cavity of 
radius r and the elastic energy required to create such a cavity. Using vf = (4-7r/3)(r — a/2) 3 , the expression for this 
free energy isPI 



3 In what follows, we leave out the factor involving the de Broglie thermal wavelength, as it cancels anyway in the final result. 
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Differentiating Eqn. ( f46| ) with respect to r yields the following equation for the equilibrium radius of the cavity: 
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This yields the following equation for r 
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{a/2 + r a )r + — 



3kT 
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and hence 



ri = 



a/2 + r , (a/2 + r ) 2 



arQ 2>kT 
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(48) 



Inserting Eqn. (^8|) in Eqn. 
lattice site 



we obtain the expression for the total free energy of an interstitial at a specific 



f I = -kTln (4V3)(r/-a/2) 3 + 8 W (r/-r ) 



(49) 



If we use the parameters for a hard-sphere solid at melting (/i w C44 = 46 [Q), wc find that the predicted concentration 
of interstitials is approximately 1 • 10~ 7 . Considering the crudeness of the approximations involved in deriving this 
result, the agreement with the corresponding numerical estimate ( xi ~ 3 • 10 -8 ) is gratifying. However, at higher 
densities, the agreement becomes worse, possibly because it is no longer justified to assume isotropic, linear, elastic 
behavior around an interstitial (see table 1) . 



In summary, we have shown that the equilibrium concentration of interstitials in hard-sphere crystals is so low 
that most experiments will not be able to detect them. We find that interstitials are quite mobile. This implies 
that interstitials that are trapped during crystal growth should be able to diffuse rapidly to the crystal surface. This 
information is good news for experimentalists trying to grow photonic bandgap materials. On the other hand, colloidal 
hard sphere crystals will have a high equilibrium concentration of vacancies. With the present accuracy of free-energy 
calculations, vacancies yield a detectable change in the free energies, but interstitials do not. 
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p 


N 


V 




gi/kT 


Xl 


h e t 


C44 


X I analytical 


11.0 


256 + 1 


0.536 


16.5 


29.9(1) 


1.5 • 10 _b 


0.087(6) 


41 


1.1 ■ 10~ b 


11.7 


256 + 1 


0.545 


17.1 


34.5(2) 


2.7 ■ 10" 8 


0.032(2) 


46 


1.3 ■ 10" 7 


11.7 


2048 + 1 


0.545 


17.1 


34.7(2) 


2.4 ■ 10~ 8 




46 


1.3 • 10~ 7 


12.0 


256 + 1 


0.548 


17.4 


36.5(2) 


5.6 ■ 10" 9 


0.079(9) 


48 


6.8 • 10~ 8 


13.0 


256 + 1 


0.559 


18.4 


44.1(3) 


7.2 • 10~ 12 


0.118(8) 


57 


3.2 • 10~ 9 



TABLE I. Simulation results for the properties of interstitials in hard sphere crystals. The values for the packing fraction 
77 and the chemical potential [i were taken from refs. 0] and js|. h e t is the fraction of interstitials found in tetrahedral holes. 
The values from the analytical estimate of section M are given as Xi t „ a i y ticai, using C44 values interpolated from fl] 
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P 

FIG. 2. xi as a function of the reduced pressures P. The drawn line corresponds to a fit of the form xi = exp(— 6. IP + 54) 
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